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FOREWORD 


This  study  represents  the  initial  effort  aimed  at  addressing  the  bend 
buckling  of  circular  cylindrical  shells  subject  to  two  self-equilibrating  end 
moments  and  zero  net  axial  loading.  Such  conditions  may  arise  in  the  flexural 
response  of  a  submerged  cylinder  excited  by  an  underwater  explosion  bubble.  The 
inertia  forces  are  neglected  here  in  view  of  the  low  frequency  nature  of  the 
response.  Both  material  and  geometrical  nonlinearities  are  included. 

This  analysis  uses  the  nonlinear  finite  element  program  ABAQUS.  Analytical 
predictions  from  ABAQUS  are  validated  by  comparison  of  results  to  experimental 
data  available  in  the  open  literature. 
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INTRODUCTION 


The  problem  of  the  deformation  of  a  straight  or  curved  shell,  subject  to 
end  bending  moments  and  possibly  internal  pressure,  has  been  addressed  in 
numerous  ar t ic les . ^ 


This  report  addresses  numerically  the  problem  of  a  straight  circular  tube 
subject  to  external  end  bending  moments  without  internal  pressure,  allowing  for 
geometric  and  material  nonlinearities.  The  nonlinear  finite  element  program 
ABAQUS^^^  is  used.  Experimental  results  from  the  open  literature^^  are 
compared  with  the  numerical  results  in  order  to  validate  ABAQUS  for  problems  of 
this  category. 


References  3,  4,  30,  and  36  are  essentially  concerned  with  the  linear 
bending  of  tubes  subject  to  end  bending  moments.  In  modern  terminology,  they 
are  referred  to  as  "geometrically  linear"  analyses.  The  first  "geometrically 
nonlinear"  analyses  are  due  to  Brazier,^  Chwalla,^^ Wood,^^  Reissner 
Keissner  and  Weinitschke , and  Weinitschke,^*^  to  mention  a  few. 


Ades,^*^  assuming  that  the  cross-sections  of  a  long  cylinder  remain 
elliptical  after  deformation,  accounted  for  geometric  and  material 
nonlinearities.  Afendik®^’®^  presented  an  approximate  analysis  incorporating 
plasticity.  References  131  and  140  allowed  for  geometrical  and  material 
nonlinearities,  while  Reference  135  extended  an  earlier  iiia  1  ■  ^  i  s  ^  to 
e lastoplastic  behavior  of  imperfect  cylinders.  References  100,  108,  115,  and 
152  are  the  only  numerical  papers  (as  applied  to  straight  cylindrical  shells)  by 
the  finite  element  method  of  which  the  author  is  aware. 


Interest  in  this  problem  stems  from  the  fact  that  during  an  underwater 
explosion,  external  vertical  forces  are  set  up  on  a  submerged  structure  and 
cause  it  to  bend,  as  if  subjected  to  end  couples.  An  analysis  could  be  of 
potential  use  to  ship  designers  when  questions  of  quantifying  ultimate 
longitudinal  strength  arise  (see  Caldwell^^).  Another  potential  problem  is 
related  to  curved  pipes  between  fixed  suppor t s . ^ ® ^ > 

TC'S^'17/^  *•  •  • 

’  ’  *  When  a  temperature  increase  occurs,  the  curved  part  is  subjected 

to  terminal  couples,  which  reduce  the  radius  of  curvature. 
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METHOD 

As  mentioned  previously,  tiie  dynamic  problem  will  be  approximated  by  a 
static  equivalent  in  view  of  the  low  frequency  of  the  motions.  The  complexity 
of  the  problem  necessitates  the  use  of  a  numerical  procedure.  The  nonlinear 
finite  element  program  ABAQUS^^^  is  used  in  this  work,. 

First,  we  establish  a  global  right-handed  coordinate  (x,y,z)  system 
(Figures  1  and  2)  with  z  the  longitudinal  axis  of  the  cylindrical  shell,  x  the 
vertical,  and  y  the  transverse  direction.  Next,  we  discretize  the  structure  by 
modeling  only  one-half  the  length  and  one-half  the  periphery,  i.e.,  we  employ  a 
quarter  model.  It  is  assumed  that  the  cylinder  is  perfectly  circular  (without 
imperfections  due  to  fabrication  and  residual  strains)  and  the  end  loads  are 
symmetrical  with  respect  to  tiie  x-z  plane,  with  the  bending  couples  lying  on  the 
global  y-axis.  The  cylindrical  surface  is  replaced  by  ABAQUS  S8R  shell  elements 
with  3  integration  points  across  the  thickness.  Geometric  and  material 
nonlinearities  are  allowed.  Perfect  plasticity,  von  Mises  isotropic  yield,  and 
an  associated  normal  flow  rule  are  used  by  ABAiJUS.  The  employed  mesh,  including 
midside  nodes,  was  25  x  25  excluding  an  auxiliary  node.  Figure  18  displays  a 
13  X  13  mesh  which  excludes  all  midnodes. 


BOUNDARY  CONDITIONS 


All  boundary  conditions  are  given  in  the  global  Cartesian  frame  of 
reference.  Along  both  generators,  AA'  (0  =  0°)  and  CC  (e  =  180°),  owing  to 
symmetry,  we  have 

V  =  <p  =  m  =0 
X  z 

Along  the  half-circle  ABC,  symmetry  implies  that  (half  length  analyzed  only) 


w  = 


<P 


X 


ip  =  0 

y 


Note  that  up  to  this  point,  the  vertical  rigid  body  motion  has  not  yet  been 
removed,  and  it  must  be  constrained  prior  to  solution.  This  will  be  done  in 
CO.  ^unction  with  the  metlii.id  of  exertin;’,  tin'  external  loading. 
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APPLIED  LOADING 


ABAQUS  has  a  particularly  attractive  feature:  a  nonlinear  multiple 
constraint  capability  (MPC).*  The  loading  on  the  structure  from  the  underwater 
explosion  is  approximated  through  an  overall  bending  moment  set  up  by  the  action 
of  longitudinal  membrane  stresses.  This  will  be  modeled  through  a  prescribed 
end  rotation  about  the  Y-axis  applied  incrementally.  For  small  deformations, 
the  vertical  coordinate  of  the  neutral  axis  from  the  center  of  the  circular 
cross-section  is  extremely  small.  Consequently,  the  sh.ift  is  approximately 
zero.  However,  for  larger  deformations  this  shift  is  substantial,  necessitating 
an  iteration  if  we  are  to  assume  that  initially  a  linear  distribution  of  forces 
produces  a  net  applied  moment.  To  avoid  the  above,  the  following  method  is  used. 

There  are  three  conditions  that  must  be  fulfilled  in  this  approach  and  are 
summarized  here  for  clarity: 

1.  Plane  sections  remain  plane  at  A'B'C'  arc  (Figure  2)  (at  the  end 
where  the  external  loading  would  have  been  applied). 

2.  There  will  be  no  rotation  of  A'B'C'  plane  (Figure  2)  about  the 
local  Y-axis. 


3.  An  end  rotation  about  the  global  y  (or  local  Z)  is  incrementally 

applied. 

Condition  3  gives  rise  to  a  distribution  of  external  longitudinal  inplane 
forces  along  the  local  X-axis,  which  causes  a  bending  moment. 

An  auxiliary  node  is  set  up  to  coincide  with  the  center  of  the  original 
undeforraed  plane.  The  terminal  couple  is  applied  at  this  section.  After 
deformation,  it  is  assumed  that  plane  sections  remain  plane.  By  St.  Venant's 
Principle  of  "elastic  equivalence  of  statically  equipollent  systems  of 
load ,  we  conclude  that  for  lengtlis  larger  than  the  cylinder  diameter, 

the  stress  distribution  away  from  the  load,  due  to  a  zero  net  axial  force  and  a 
terminal  bending  moment,  does  not  depend  on  the  traction  distribution,  except 
perhaps  locally  in  the  neighborhood  of  the  point  of  application  of  the  load. 
This  condition  can  be  fulfilled  if  vector  £3  (Figure  2)  of  the  local  frame  of 
reference,  located  on  the  deformed  plane,  is  orthogonal  to  any  vector  on  that 
plane. 

The  position  vectors  before  deformation  of  nodes  on  the  periphery  of  the 
shell,  where  the  end  couple  is  applied,  are  £(,,  and  for  the  auxiliary 

node.  After  displacements  u  of  the  peripheral  nodes  and  of  the  auxiliary 
node,  we  obtain 


r  =  r  +  u 
—  — o  — 


(1) 


b  D  b 

r  =  r  +  u 
—  — o  — 


(2) 


*MPC  =  Multiple  Point  Constraints 
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br  bif 

r-r  ■Ir^-r^l+lu-u] 


Therefore , 


*  I3  =  0 

where,  in  component  notation, 


r 


((x  -  Xfa),  (y  -  iz  -  z^)] 


and  all  quantitites  are  given  with  respect  to  the  global  undeformed  frame, 
if  the  original  local  x-axis  (along  vector  had  an  initial  inclination 
(Pq  with  respect  to  the  global  y-axis,  after  deformation  the  angle  would  be 
•^o  Therefore,  after  deformation  the  direction  cosines  of  the 

local  vectors  i  £2  *  ^  with  respect  to  the  global  system  would  be 


=  [cos((p^  +  (Py^),  0,  -sin((p^  + 


a  =  [0,  1,  0] 


a»  =  [sin((p  +  ;p  ^),  0,  cos(q)  +  <f  '’)] 

— j  o  y  o  y 


Therefore , 


(x 


X,  )sin(<p  + 
b  o 


) 


b 

+  (z  -  z  )cos((p  +  (p  )  =  0 
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This  can  be  further  expanded  by  means  of  the  trigonometric  identities 


cos(<B  +  m  ) 

o  y 


b  .  .  b 

COSO  cos<p  -  sino  sintp 
o  y  o  y 


(10) 


sin(<B  +0  ) 

o  y 


b  ,  b 

Sin<p  COS()>  *  C08<P  Blivp 
o  y  o  y 


(11) 


Since  our  solution  process  is  incremental  and  our  constraints  are  nonlinear,  our 
boundary  conditions  must  be  cast  in  incremental  form.  Consequently,  we  note 
that  <Po  is  a  constant  angle  (which  for  straight  tubes  is  zero).  The 
incremental  form  of  the  left-hand  side  of  Equation  (3)  can  be  obtained  from  the 
difference  of  the  incremental  forms  of  Equations  (1)  and  (2),  since 


r+4r=r  +u+Au 
—  —  — o  —  — 


(r  +  Ar)  - 

+  (Au  - 
In  view  of  Equation 
Ar  -  Ar  = 


b  b  b 
r  +  u  +  Au 
— o  —  — 


(£^  +  Ar^)  =  (r^ 


Au*") 

(3). 

Au  -  Au' 


r  *")  +  (u  -  u^*) 
— o  —  _ 


(12) 

(13) 


(14) 


(15) 


The  incremental  form  of  the  righthand  side  of  Equation  (9)  with  respect  to 
increments  of  the  displacements  U2,  u^,  ?y^,  Ug*’,  and  Uj^^,  by  means  of 
Equation  (15)  becomes 


(Au  -  Au  ^)[cos<p  cos9  ^  -  sin9rvSin9 
z  z  o  y  y 


+  (Au  -  Au  ^)[sin<p  C089  ^  +  cos®  sin9 
XX  o  y  o  y 


(16) 


+  ( z 


z,)l-cos<p  sin<p 
b  o  y 


b 

sin®  cos® 
o  y 


lA® 


b 

y 


+  (x  -  X,  )  [-sin®  sin® 


b,  b 

+  cos®  cos®  JA®  =  0 
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or 


(Au  -  Au  )cos((p  + 
z  z  o 


<P  *") 

y 


+  (Au  -  Au  )sin((p  + 
X  X  o 


<P  ) 

y 


+  Aip*^ 

y 


z.  )sin(<p  +  (p 
b  o  y 


(x 


X,  )cos((p  + 
b  o 


0 


(1- 


Next,  we  note  the  end  plane  does  not  rotate  about  the  local  Y-axis,  i.e.,  thei 
will  be  no  component  of  rotation  about  vector  A  small  rotation  about  the 

local  Y-axis  can  be  obtained  by  a  rotation  vector  with  respect  to  the  global 
system 


-  ^‘*'x’  ‘*’y’  '’’z^ 

If  this  vector  is  normal  to  vector  a^ ,  then  it  has  no  rotational  component 
about  ,  i.e., 


£  •  a  =  (<p  ,  <p  ,  (p  ]  •  [cos((p  +  (p  ),  0 

— i  X  y  z  y  y 


-sinC'cp  +  <p  )]  =  0 
o  y 


or 


<p  co8(<p  +  ip  ^)  -  <p  sin((p  +  9  ^)  =  0 
X  ’^o  y  z  o  y 


The  incremental  form  of  Equation  (20)  becomes 


A<p  cosCffl  +  (p  -  A<p  8in(<p  +  <p 
X  o  y  z  o  y 


-  A(p  ^  |(p  sinCm  +  <p  ^)  +  (p  cosCip  +  9  ^))  =  0 


(1< 


(2C 


(2 
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Finally,  the  essential  boundary  condition  is  that  the  rotation  of  all 
nodes  on  the  periphery  must  equal  the  rotation  about  the  local  Z-axis  ( a? 
vector  or  global  y-axis)  of  the  auxiliary  node,  i.e.. 


and  in  incremental  form 


Aq)  =  A<p  (23: 

y  y 

The  vertical  (global  x-axis)  rigid  body  motion  must  be  removed.  This  can  be 
accomplished  by  coupling  the  average  X-displacement  of  the  shell  nodes  on  the 
end  section  to  the  vertical  displacement  of  the  auxiliary  node.  This,  in  turn 
is  set  to  zero.  For  element  i,  S8R  element  displacements  across  an  edge  are 
quadratic,  i.e.  (with  0  ^  5  £  1)  the  vertical  displacement  u(5)  is  given  in 
terms  of  its  nodal  values  (uj^,  uv^+i,  uk+2^ 


niO  =  [(2^  -  1)(5  -  1),  4C(1  -  C),  t(2C  -  1)] 


\ 


k+1 


1-^2  J 


(24: 


For  uniform  finite  element  grid,  where  s  denotes  arc  length,  s£  arc  length  of 
element  i,  and  L  total  arc  lengthi 


ds.  =  4.dC  =  r  d?  (26: 

II  N 

At  the  element  level  (element  i),  with  nodal  vertical  displacements  (uy^, 

“k+l>  “k+Z^ 


n 


NSWC  TR  84-380 


8“8  . 


f  ‘  T  r 

J  u(t)d8.  »  /  u(5)l.d5  =  i.  /  U 

1  i-'a 


(Ode 


8“0 


\  ■ 

1 . 

*  ^  [1,  4,  1) 

'^k  +  1 

6N  ^ 

\  +  l 

\+2 

\+2 

Since  there  are  N  elements  (with  2  corner  and  1  midside  nodes,  i.e.  3  nodes 
total),  there  will  be  2N+1  nodes  and  associated  nodal  values  of  vertical 
displacements.  Therefore,  by  Equation  (24) 


s=L 


i=N  5  =  1 


s=0 


J"  u(s)ds  =  ^  ^  J'  u(5)t. 

A  »-- />  ^ 


de 


i=l 


e=o 


4,  2, 


4, 


{lx(2N+l) } 


2,  4,  11 


3+1 


S+2N-1 


S+2N 
{(2N+1)  X  1) 


Consequently,  the  boundary  condition 


L 


s=L 


u (s )ds 


b 

u 

X 


can  be  recast,  in  view  of  Equation  (28),  to 
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6N 


[1.  4. 


2.  4. 


2,  4.  1) 


8+1 


8+2N-1 


b 

u 

X 


{1  X  (2N+1)} 


“8+2N 


{(2N+1)  X  1} 


(30) 


for  a  uniform  grid,  where  the  nodal  vector  of  Equation  (30)  denotes  vertical 
displacements  of  shell  nodes  at  the  end  plane.  The  above  constraints  can  be 
incorporated  through  the  nonlinear  MFC  capability  of  ABAQUS.^52  Furthermore, 
from  Equations  (3)  and  (9),  Equation  (9)  can  be  rewritten  as 


u 

X 


b 

X 

o 


]sin(<p^ 


+ 


) 


(31) 


+  [z  +  u 
o  z 


u  *^]cos(9  +  9  ^) 

2  o  y 


0 


Following  Reference  132,  we  eliminate  u^  from  Equation  (31)  first  (i.e.  at  the 
shell  node) 


u 

z 


\  /  ^ 
z)-vx  -  X  +u 

O  O  O  X 


u  '^)tan(9  +  9  ^) 

X  o  y 


Solving  for  9^  from  Equation  (20) 


(32) 


9=9  tan(9  +  9  ^)  (33) 

X  z  o  y 


and  9y  from  Equation  (22) 


9 


y 


(34) 


Furthermore,  we  eliminate  the  vertical  translation  u^.  This  follows  from 
constraints  involving  the  elimination  of  u^.  All  previous  steps  are  explained 
in  Reference  152  or  can  be  found  in  the  FORTRAN  listing  for  MFC  constraints. 
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APPLICATIONS  AND  RESULTS 


Four  unstiffened  circular  cylinders  have  been  analyzed  using  the  finite 
element  program  ABAQUS^^^:  models  lOA,  16A,  20A  of  Reference  93,  and  20AI, 
which  is  an  imperfect  version  of  model  20A.  Note  that  they  are  comparatively 
thick  to  avoid  premature  collapse  and  to  exhibit  plasticity  effects  as  the 
applied  external  moment  increases  past  the  yield  moment.  The  above  models  fall 
in  the  range  of  "long  cylinders."  Their  perfect  discrete  analogues  would  fail 
by  ovalization  or  the  "Brazier  effect,"  or  by  plastic  deformation.  The 
experimentally  imperfect  models  failed  by  buckling  "failure"^^  as  displayed  on 
the  moment-curvature  curves. 

Table  1  gives  the  geometrical  and  material  properties  of  these  models  and 
the  stress-strain  curves  can  be  obtained  from  Reference  93.  In  Reference  93, 
the  external  loading  was  applied  by  the  use  of  a  so-called  "shear  span"  prior  to 
the  test  section  "bending  span."  To  create  the  external  moment  in  our  research, 
we  employed  an  additional  span  beyond  the  bending  span,  equal  to  the  "shear" 
span,  and  then  applied  a  fixed  angle  of  rotation.  Table  2  contains  some 
parameters  used  in  reducing  the  stresses,  moments,  curvatures,  and  forces  in 
nondimens ional  form,  and  they  can  be  used  to  determine  relative  magnitude  for 
critical  quantities  such  as  yield  moment,  etc. 

Mean  curvature  k  is  defined  as  the  ratio  of  the  sum  of  the  absolute  values 
of  direct  longitudinal  strains  at  0*  (top  of  cylinder,  i.e.  tension  side)  and 
180“  (bottom  of  cylinder,  i.e.  compression  aide)  divided  by  the  undeformed 
diameter  of  the  shell.  Figures  3  through  5  are  the  relevant  M  vs.  k  curves  for 
models  lOA,  16A,  and  20A,  together  with  the  corresponding  experimental  results 
of  Reference  93.  Agreement  between  experimental  (dots)  and  computed  results  is 
good  for  all  three  cylinders.  Imperfection  sensitivity  and  the  analysis  of 
medium  and  short  cylinders,  where  short-wave  length  buckling  may  control  the 
collapse  mechanism,  will  be  addressed  elsewhere. 

Figures  6  through  8  display  the  work  done  versus  the  moment  parameter  u  = 
Mext/^0>  where  Mg  =  4R^  hOy,  R  =  mean  radius  of  cylinder,  E  =  Young's  Modulus, 
Mrxt  “  external  bending  moment,  h  =  shell  thickness,  and  Oy  =  material  yield 
stress . 

We  define  two  parameters,  longitudinal  stress/stress  at  bifurcation  and 
hoop  stress/stress  at  bifurcation,  where  stress  at  bifurcation  is 


=  E  h 

(l-v2)  ^ 

Figures  9  through  11  represent  plots  of  the  longitudinal  inplane  stress 
parameter  versus  angular  position.  Additional  local  bending  stresses  across  the 
shell  thickness  are  not  addressed.  Note,  however,  that  the  longitudinal 
membrane  stress  distribution  in  Figures  9  through  11  does  not  follow  simple  beam 
theory.  The  stress  distribution  from  the  0“  and  180“  points  (top  and  bottom  of 
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EXTERNAL  WORK  -  SENDING  MOMENT/M.jlthvi  MOMENT 


0  20  40  60  80  100  120  140  160 

WORK  DONE  (KIP  IN) 

FIGURE  6.  M/Myijliy,  WORK  DONE  (KiP  IN)  FOR  MODEL  10A 


FIGURE  9.  LONGITUDINAL  STRESS  (STEP  17)  VERSUS  ANGULAR  POSITION  FOR  MODEL  10A 


FIGURE  10.  LONGITUDINAL  STRESS  (STEP  31)  VERSUS  ANGULAR  POSITION  FOR  MODEL  16A 


FIGURE  11.  LONGITUDINAL  STRESS  (STEP  20)  VERSUS  ANGULAR  POSITION  FOR  MODEL  20A 


NSWC  TR  84-38C 


the  iial  t-sec  t  ion )  is  reduced  (.tension  zone)  or  increased  (compression  zone) 
slightly  betore  it  as.suines  a  linear  torra. 

Figures  12  through  14  display  the  variation  ot  the  hoop  stress  parameter 
with  angular  location.  The  distribution  ot  this  stress  over  the  lia I f -sec t ion 
clearly  indicates  maximum  compressive  hoop  stresses  around  the  90“  location, 
with  corresponding  maximum  tensile  hoop  stresses  at  0°  and  180°,  respectively. 
These  stresses,  however,  are  smaller  than  the  longitudinal  membrane  stresses  by 
an  order  of  magnitude.  All  stresses  shown  correspond  to  the  final  point  of  the 
moment-curvature  plots. 

Close  examination  of  Figures  15  through  17  reveals  that  the  M/M(;g  vs.  R/kcR 
curves  have  a  slope  at  the  origin  of  approximately  1.0.  This  agrees  fairly  well 
witn  the  initial  slope  predicted  by  Reissner's  nonlinear  theory^^  as  well  as 
Von  Karman's  linear  analysis.^  In  the  present  notation,  Reissner's  relation¬ 
ship  between  M/M^r  and  R/RcR  written  as 


=  k/k^j^[1.0  -  0.5  -  (1/6)  (k/k^^)^  .  .  .  .  J  (35) 


with  an  obvious  slope  of  1.0  at  the  origin,  maximum  value  of  M/M(jR  =  0.5011 
at  k/kcR  =  0.7  1954.  The  maximum  value  of  0.5011  is  much  higher  than  both 
experimental  and  numerical  results,  showing  the  effects  of  plasticity  which 
Reissner's  theory  does  not  account  for.  Note  that  plots  of  p  =  (MgxT^^O^ 
versus  work  done  (Figures  6  through  8)  clearly  display  that  ultimate  moment  is 
reached  at  u  =  1  as  predicted  by  beam  plastic  theory.  Actually,  it  exceeds 
P  =  1  by  about  10  percent. 

Figure  18  displays  the  deformed  and  undeformed  profiles  of  cylindrical 
model  20A.  Figure  19  represents  the  corresponding  deformed  and  undeformed 
profiles  of  the  imperfect  version  of  model  20A,  which  from  now  on  will  be 
referred  to  as  model  20A1.  Figure  20  is  an  enlargement  of  the  cross-sections  at 
midlength  before  and  after  deformation  of  model  20AI.  Figure  21  shows  the 
moment-curvature  plots  of  model  20AI  compared  with  the  experimental  data 
(circled  points).  Similarly,  Figures  22  through  25  show  results  for  model  20AI 
of  the  following  pairs  of  variables: 


p  =  MlXt/^O  ^/^ULTIMATE^  ~  Work  Done  (kip-in.) 


LONG 

o  /o  ~ 

INPLANE  CK 


Angular 


Position 


j  hOUPyjj  Angular  Position  and  M/Mcr  ~  k/krg 


CAL 
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pijpt. 

im.  ricTon  •  m.V'H 
snio  L1MC9  '  OlSPL^D  NCIH 
OMHCD  IIHC9  (MKtMH.  MCSH 


1  ^srcf  1  IMdCnCRT  59 


FIGURE  18.  CYLINDRICAL  SHELL  (MODEL  20A)  SUBJECT  TO  END  BENDING  MOMENT  VIEWED  FROM 
100  MOO  ",500  "  AT  STEP  3,  INCREMENT  55 
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■'A. 

"n. 

■■■CfLINDIR  WITH  IMFcRFt!  ’’li.Tj-.  ir.ifll  TO  i:Nj  TOUFLE 

t  ]t«C«CNClT  2 


figure  19.  IMPERFECT  CYLINDRICAL  SHELL  {VARIANT  OF  MODEL  20A)  SUBJECT  TO  END  BENDING 
MOMENT  VIEWED  FROM  100  MOO  '. 500'  AT  STEP  6  INCREMENT  2 
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FIGURE  20, 


UNDEFORMED  AND  DEFORMED  HALF  CROSS  SECTION  OF  IMPERFECT  VERSION  OF 
MODEL  20A  AT  MIDLENGTH 
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BENDING  MOMENT  -  MEAN  CURVATURE 


FIGURE  21.  MOMENT-CURVATURE  PLOTS  FOR  IMPERFECT  VERSION  OF  MODEL  20A 


EXTERNAL  WORK  -  BENDING  MOMENT/ULTIM  MOMENT 
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FIGURE  22.  M/MULTIM  -  WORK  DONE  (KIP-IN.)  FOR  IMPERFECT  VERSION  OF  MODEL  20A 
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FIGURE  23.  LONGITUDINAL  STRESS  (STEP  6.  INCREMENT  2,  WHICH  CORRESPONDS  TO  FINAL 

PLOTTED  POINT  ON  FIGURE  21)  VERSUS  ANGULAR  POSITION  FOR  IMPERFECT  VERSION 
OF  MODEL  20A 


FIGURE  24.  HOOP  STRESS  DISTRIBUTION  (STEP  6,  INCREMENT  2,  WHICH  CORRESPONDS  TO  FINAL 
PLOTTED  POINT  ON  FIGURE  21)  VERSUS  ANGULAR  POSITION  FOR  IMPERFECT  VERSION 
OF  MODEL  20A 
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At  this  point,  note  that  the  relevant  figures  corresponding  to  the  perfect  model 
20A  are  8,  11,  14,  and  17,  respectively.  The  imperfection  in  the  radial 
displacement  of  the  original  surface  from  the  mean  radius  R  of  imperfect  model 
20AI  varied  according  to  the  formula 


h  sin  cos (10  0) 

Imperfect  model  20AI  was  generated  with  one-half  wave  axially  and  10  waves 
peripherally . 

Figure  26  presents  superimposed  moment-curvature  curves  for  both  models  20A 
and  20AI.  Notice  that  since  they  are  comparatively  thick  and  fail  by 
plastification,  unlike  ovalization  or  bifurcation  failure,  they  are  not 
imperfection  sensitive.  The  response  of  the  imperfect  model  20AI  is  very 
similar  to  the  perfect  one. 

Table  3  summarizes  digitized  results  pertaining  to  bending  moment,  angle  of 
rotation  at  the  end,  where  the  external  load  is  applied,  as  well  as  strain 
energy,  work  done,  and  plastic  dissipation  for  model  20A.  Note  that  because  of 
an  existing  error  in  ABAQUS,  concerning  how  energies  are  computed,  strain  energy 
and  plastic  dissipation  do  not  agree  exactly  with  the  work  done  (0.120489  + 
0.290966  =  0,411455  compared  with  0.417864).  In  fact,  this  difference  decreases 
as  we  march  along  the  load-deformation  curve.  Table  4  gives  the  stress 
distribution  for  load  step  20  of  model  20A.  Tables  5  and  6  give  the 
corresponding  information  for  model  20AI,  and  Tables  7  and  8  give  the  actual 
points  plotted  in  Figure  26  for  both  models  20A  and  20A1 . 

Finally,  a  fine  and  extremely  important  point  pertaining  to  the  modeling 
issue  must  be  addressed.  The  experimental  set  up  involved  the  analysis  of  a 
three  span  beam  (in  all  cases)  subject  to  vertical  self-equilibrating  shear 
loads,  to  simulate  overall  bending  over  the  two  end  spans.  The  middle  span  had 
no  loads  applied.  In  this  analysis,  one  end  span  and  half  of  the  center  span 
(symmetry)  were  modeled.  A  rotation  was  enforced  through  the  "auxiliary  node" 
concept  and  the  MFC  constraints.  Subsequent  computations  without  the 
"additional"  end  span  gave  a  slightly  different  response. 

Consequently,  what  constitutes  an  adequate  additional  span  to  be  included 
for  proper  response  has  not  been  determined  but  is  of  major  importance  in  this 
analysis.  In  addition,  to  complete  studies  on  the  bend-buckling  modeling  of 
cylindrical  shells,  "short"  and  "medium"  length  tubes  with  or  without  ring 
stiffeners  must  be  addressed. 

In  closing  this  discussion  it  must  be  stressed  that,  by  neglecting  inertia 
effects,  a  reasonably  successful  first  approximation  of  the  problem  of  a  long 
straight  circular  cylinder  shell  subject  to  end  couples  has  been  developed. 
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FIGURE  26.  bending  MOMENT  DISTRIBUTION  VERSUS  MEAN  CURVATURE  FOR  BOTH  PERFECT 
AND  IMPERFECT  VERSIONS  OF  MODEL  20A 
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TABLE  8.  MOMENT-CURVATURE  RESULTS  FOR  IMPERFECT  VERSION  OF  MODEL  20A  (MODEL  20At) 


CURVATURE  K 
(1/FT) 

U.  UUS'PSOUUbrS© 
0.  3379i.i;95E-03 
0 . 15754i41E-0£ 
0,  ?£l£vi5o0E~3E 
0 . 3150431 3E-0E 
2  .  5  3 1 9  9 1  6  8  E  “  5 
0  5  4-  5  9  3  5  3  8  E  “■  0  E 
S . 36537399E-33 
3 . 333S33S9E-03 
0  3  3  9  4  3  4  7  0  E  ~  0  3 
0. 5yu3v559E-03 
0 . 41S3351 lE-03 
0 , 439170i3E-03 
0. 4457094SE-03 
0.  4iv9?4S27E-03 
0 . 4737371 lE-02 
0 .  495  33E5  9E-02 
3. 53174359E-0a 


MOMENT  M 
(KIP-FT) 

0  0C000000E+00 
0 .  5  5431 1 21E  +  05 
0  .  197649i.0E  +  0  3 
0 . 337950113+03 
0  .  3  9  9  4  5  3  9  5  E  +  0  3 
0. 405~4533£+03 
0. 41 093304E+03 
0. 4155 i355E+03 
3 . 4195935 lE+03 
0 . 45093170E+03 
0. 45509 50 EE+03 
0. 4545 151 lE+03 
0. 45539 ITSEl 05 
0. 459i6354E+03 
0. 431ii319E+03 
0. 4255t-930E  +  O3 
0. 43395i39E+03 
0. 4344110iE+03 
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SUMMARY 


A  modeling  strategy  is  established  to  obtain  the  moment-curvature  relation 
as  well  as  the  relation  between  the  work  done  and  the  applied  moment  for 
circular  cylindrical  shells.  This  is  achieved  by  using  the  nonlinear  finite 
element  program  ABAQUS  in  conjunction  with  preprocessing  and  postprocessing 
computer  programs.  The  results  compare  favorably  with  experimental  curves 
reported  in  the  open  literature.  Such  analysis  is  of  potential  use  in 
predicting  critical  bending  moments,  ultimate  moment  for  ship  hulls,  pipe  bends 
in  nuclear  reactors,  submarine  pipelines,  etc. 
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E 

FO 

h 

I 

k 


L 


^br 

Mcr 

M,1%xT 


Mq 

N 


=  Local  vectors  along  the  local  Y,  Z,  and  X  axes.  Local  z-axis 
is  parallel  to  global  y.  Initially,  local  y  is  not  parallel, 
but  is  coplanar  with  global  x  and  in  the  same  direction. 
Initially,  local  X  is  coplanar  with  global  z  and  points  in  the 
same  direction  (Figure  2). 

=  Young's  Modulus 

=  Axial  force  based  on  yield  stress 

=  Critical  curvature  at  bifurcation  (see  Table  2) 

=  Shell  thickness 

=  Moment  of  inertia  of  undeformed  section  (Table  2) 

=  Subscript  in  nodal  values  uj^,  ur+i,  uj^+2  (Eq.  (24)). 

Mean  curvature  defined  as  the  sum  of  absolute  value  of  direct 
axial  strains  at  the  top  and  bottom  of  the  cylinder  divided  by 
the  undeformed  diameter  of  the  shell. 

=  Total  length  of  cylindrical  shell  (Figure  1).  In  Equation 

(25),  L  is  total  arc  length  of  end  section  ABC  (Figure  2)  prior 
to  deformation. 

=  Length  of  side  of  element  i  on  arc  where  end  rotation  is 
applied  (undeformed  side  ABC  on  Figure  2) 

=  Critical  bending  moment  by  Brazier  (Table  2) 

=  Critical  bending  moment  of  bifurcation  (Table  2) 

=  External  bending  moment  calculated  in  the  form  of  a  reaction  in 
this  analysis 

=  Plastic  moment  or  ultimate  moment  based  on  yield  stress 
(Table  2) 

=  Total  number  of  equal  elements  on  undeformed  arc  ABC  (Figure  2) 
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£  =  Final  position  vector  on  arc  A'B'o'  (Figure  2,  after 

de  formation) 

Tq  =  Initial  position  vector  o*"  any  point  on  periphery  of  shell  (arc 

ABC  on  Figure  2)  at  which  the  end  rotation  will  be  applied 

=  Initial  position  vector  of  auxiliary  node 


R 


=  Mean  shell  radius 


u  =  In  the  "APPLIED  LOADING"  section,  u  is  the  vector  displacement 

after  deformation;  it  is  not  to  be  confused  with  u 

=  In  the  context  of  analysis  given  in  "APPLIED  LOADING,"  vector 
displacement  of  auxiliary  node  after  deformation;  it  is  not  to 
be  confused  with  u 


Ujf  ,  Uy  ,  u^ 


=  Vertical,  transverse  and  longitudinal  translations  (along 
global  X,  Y,  Z  axes).  Also  referred  to  as  u,  v,  and  w. 


(u) 


=  Vertical  displacement  at  any  point  C  on  an  element  side 


Vl>V2 


=  Nodal  values  of  displacement  u(^)  at  end  node  k,  midside  node 
k+1,  end  node  k+2. 


*^8 

u 

V 

V 

0 

V 


=  Nodal  value  at  s  (point  on  arc  where  end  rotation  is  being 
applied) 

=  Translation  along  global  x-axis 
=  Translation  along  global  y-axis 
=  Translation  along  global  z-axis 

=  Angle  from  vertical  global  x-axis  up  to  end  of  arc  to  be 
analyzed 

=  Poisson's  ratio 


K  =  Local  normalized  variable  in  interval  (0,  l)  employed  in 

describing  vertical  displacement  distribution  u(c)* 

=  Critical  stress  at  bifurcation  (Table  2) 

oy  =  Yield  stress  of  material 

<pQ  =  Initial  inclination  of  cylindrical  shell  about  y-axis.  For 

straight  tubes  (present  case)  fo  “  ® 

Vx  =  Rotation  about  global  x-axis 
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FOREWORD 


This  study  represents  the  initial  effort  aimed  at  addressing  the  bend 
buckling  of  circular  cylindrical  shells  subject  to  two  self-equilibrating  end 
moments  and  zero  net  axial  loading.  Such  conditions  may  arise  in  the  flexural 
response  of  a  submerged  cylinder  excited  by  an  underwater  explosion  bubble.  The 
inertia  forces  are  neglected  here  in  view  of  the  low  frequency  nature  of  the 
response.  Both  material  and  geometrical  nonlinearities  are  included. 

This  analysis  uses  the  nonlinear  finite  element  program  ABAQUS.  Analytical 
predictions  from  ABAQUS  are  validated  by  comparison  of  results  to  experimental 
data  available  in  the  open  literature. 
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INTRODUCTION 


The  problem  of  the  deformation  of  a  straight  or  curved  shell,  subject  to 
end  bending  moments  and  possibly  internal  pressure,  has  been  addressed  in 
numerous  ar t ic les . ^ 

This  report  addresses  numerically  the  problem  of  a  straight  circular  tube 
subject  to  external  end  bending  moments  without  internal  pressure,  allowing  for 
geometric  and  material  nonlinearities.  The  nonlinear  finite  element  program 
ABAQUS^^^  is  used.  Experimental  results  from  the  open  literature^^  are 
compared  with  the  numerical  results  in  order  to  validate  ABAQUS  for  problems  of 
tliis  category. 


References  3,  4,  30,  and  36  are  essentially  concerned  with  the  linear 
bending  of  tubes  subject  to  end  bending  moments.  In  modern  terminology,  they 
are  referred  to  as  "geometrically  linear"  analyses.  The  first  "geometrically 
nonlinear"  analyses  are  due  to  Brazier,^  Chwalla  Wood,^^  Reissner 

Reissner  and  Weinitschke , and  Weinitschke to  mention  a  few. 


Ades,^®  assuming  that  the  cross-sections  of  a  long  cylinder  remain 
elliptical  after  deformation,  accounted  for  geometric  and  material 
nonlinearities.  Afendik®^>°^  presented  an  approximate  analysis  incorporating 
plasticity.  References  131  and  140  allowed  for  geometrical  and  material 
nonlinearities,  while  Reference  135  extended  an  earlier  analvsis^^^  to 
e lastoplas t ic  behavior  of  imperfect  cylinders.  References  100,  108,  115,  and 
152  are  the  only  numerical  papers  (as  applied  to  straight  cylindrical  shells)  by 
the  finite  element  method  of  which  the  author  is  aware. 


Interest  in  this  problem  stems  from  the  fact  that  during  an  underwater 
explosion,  external  vertical  forces  are  set  up  on  a  submerged  structure  and 
cause  it  to  bend,  as  if  subjected  to  end  couples.  An  analysis  could  be  of 
potential  use  to  ship  designers  when  questions  of  quantifying  ultimate 
longitudinal  strength  arise  (see  Caldwell^^).  Another  potential  problem  is 
related  to  curved  pipes  between  fixed  suppor  t  s  .  ^  ®  ^  *  ^*^  *^^  * 

When  a  temperature  increase  occurs,  the  curved  part  is  subjected 
to  terminal  couples,  which  reduce  the  radius  of  curvature. 
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METHOD 

As  mentioned  previously,  the  dynamic  problem  will  be  approximated  by  a 
static  equivalent  in  view  of  the  low  frequency  of  the  motions.  The  complexity 
of  the  problem  necessitates  the  use  of  a  numerical  procedure.  The  nonlinear 
finite  element  program  ABAQUS^^^  is  used  in  this  work. 

First,  we  establish  a  global  right-handed  coordinate  (x,y,z)  system 
(Figures  1  and  2)  with  z  the  longitudinal  axis  of  the  cylindrical  shell,  x  the 
vertical,  and  y  the  transverse  direction.  Next,  we  discretize  the  structure  by 
modeling  only  one-half  the  length  and  one-half  the  periphery,  i.e.,  we  employ  a 
quarter  model.  It  is  assumed  that  the  cylinder  is  perfectly  circular  (without 
imperfections  due  to  fabrication  and  residual  strains)  and  the  end  loads  are 
symmetrical  with  respect  to  t'ne  x-z  plane,  with  the  bending  couples  lying  on  the 
global  y-axis.  Tlie  cylindrical  surface  is  replaced  by  ABAQUS  S8R  shell  elements 
with  3  integration  points  across  the  thickness.  Geometric  and  material 
nonlinearities  are  allowed.  Perfect  plasticity,  von  Mises  isotropic  yield,  and 
an  associated  normal  flow  rule  are  used  by  ABA(^US.  The  employed  mesh,  including 
midside  nodes,  was  25  x  25  excluding  an  auxiliary  node.  Figure  18  displays  a 
13  x  13  mesh  which  excludes  all  raidnodes. 


BOUNDARY  CONDITIONS 


All  boundary  conditions  are  given  in  the  global  Cartesian  frame  of 
reference.  Along  both  generators,  AA'  (9  =  0“)  and  CC  (6  =  180®),  owing  to 
symmetry,  we  have 


V  =  m  =  m 

X  Z 


0 


Along  the  half-circle  ABC,  symmetry  implies  that  (half  length  analyzed  only) 


w  =  o  =  ip  =0 
X  y 

Note  that  up  to  this  point,  the  vertical  rigid  body  motion  has  not  yet  been 
removed,  and  it  must  be  constrained  prior  to  solution.  This  will  be  done  in 
CO.  _unction  with  the  metliod  of  exerting,  the  external  loading. 
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ABAQUS  has  a  particularly  attractive  feature;  a  nonlinear  multiple 
constraint  capability  (MPC).*  The  loading  on  the  structure  from  the  underwater 
explosion  is  approximated  through  an  overall  bending  moment  set  up  by  the  action 
of  longitudinal  membrane  stresses.  This  will  be  modeled  through  a  prescribed 
end  rotation  about  the  Y-axis  applied  incrementally.  For  small  deformations, 
the  vertical  coordinate  of  the  neutral  axis  from  the  center  of  the  circular 
cross-section  is  extremely  small.  Consequently,  the  sh.ift  is  approximately 
zero.  However,  for  larger  deformations  this  shift  is  substantial,  necessitating 
an  iteration  if  we  are  to  assume  that  initially  a  linear  distribution  of  forces 
produces  a  net  applied  moment.  To  avoid  the  above,  the  following  method  is  used. 

There  are  three  conditions  that  must  be  fulfilled  in  this  approach  and  are 
summarized  here  for  clarity: 

1.  Plane  sections  remain  plane  at  A'B'C'  arc  (Figure  2)  (at  the  end 
where  the  external  loading  would  have  been  applied). 

2.  There  will  be  no  rotation  of  A'B'C'  plane  (Figure  2)  about  the 
local  Y-axis. 

3.  An  end  rotation  about  the  global  y  (or  local  /.)  is  incrementally 

applied . 

Condition  3  gives  rise  to  a  distribution  of  external  longitudinal  inplane 
forces  along  the  local  X-axis,  which  causes  a  bending  moment. 

An  auxiliary  node  is  set  up  to  coincide  with  the  center  of  the  original 
undeforraed  plane.  The  terminal  couple  is  applied  at  this  section.  After 
deformation,  it  is  assumed  that  plane  sections  remain  plane.  By  St.  Venant's 
Principle  of  "elastic  equivalence  of  statically  equipollent  systems  of 
load , we  conclude  that  for  lengths  larger  than  the  cylinder  diameter, 
the  stress  distribution  away  from  the  load,  due  to  a  zero  net  axial  force  and  ii 
terminal  bending  moment,  does  not  depend  on  the  traction  distribution,  except 
perhaps  locally  in  the  neighborhood  of  the  point  of  application  of  the  load. 

This  condition  can  be  fulfilled  if  vector  ^  (Figure  2)  of  the  local  frame  of 
reference,  located  on  tlie  deformed  plane,  is  orthogonal  to  any  vector  on  that 
plane . 


The  position  vectors  before  deformation  of  nodes  on 
shell,  where  the  end  couple  is  applied,  are  £q,  and  f 

node.  After  displacements  u  of  the  peripheral  nodes  and 
node,  we  obtain 


the  periphery  of  the 
or  the  auxiliary 
u^  of  the  auxiliary 


r  =  r  +  u 
—  — o  — 


(1) 


b  b  b 
r  =  r  +  u 
—  — o  — 


(2) 


*MPC  =  Multiple  Point  Constraints 


V  S 
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8  r  8l,r  8^ 

r-r  “[r^-r^J+lu-u] 

Therefore , 

[r  -  r^]  *  =  0 

where,  in  component  notation, 

r  -  r^  =  [  (x  -  x^) ,  (y  -  y^^),  (z  -  z^)] 

and  all  quantitites  are  given  with  respect  to  the  global  undeformed  frame, 
if  the  original  local  x-axia  (along  vector  had  an  initial  inclination 
Vq  with  respect  to  the  global  y-axis,  after  deformation  the  angle  would  be 
1*0  <Py^*  Therefore,  after  deformation  the  direction  cosines  of  the 

local  vectors  £2,  and  a^  with  respect  to  the  global  system  would  be 
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This  can  be  further  expanded  by  means  of  the  trigonometric  identities 


cosCo  + 
o 


b  .  .  b 

co8q)^co8<p^  -  sintp^sinq) 


(10) 


sin((p  + 
o 


b  b 

sino  cos<fi  *  COSO  sino 
o  y  o  y 


(11) 


Since  our  solution  process  is  incremental  and  our  constraints  are  nonlinear,  our 
boundary  conditions  must  be  cast  in  incremental  form.  Consequently,  we  note 
that  <Po  is  a  constant  angle  (which  for  straight  tubes  is  zero).  The 
incremental  form  of  the  left-hand  side  of  Equation  (3)  can  be  obtained  from  the 
difference  of  the  incremental  forms  of  Equations  (1)  and  (2),  since 


r^  +  Ar^  =  li 


b  b  b  b  b 
r  +Ar  =  r  +u  +Au 
—  —  — o  —  — 


(£  +  A£)  -  (£^  +  A£^)  =  (£^  -  £^^)  +  (u  -  u^) 


+  (Au  -  Au^) 


In  view  of  Equation  (3), 


b  I 

Ar  -  Ar  =  Au  -  Au 


(12) 

(13) 


(14) 


(15) 


The  incremental  form  of  the  righthand  side  of  Equation  (9)  with  respect  to 
increments  of  the  displacements  Ug,  u^,  "fy^*  '^x^*  “eans  of 

Equation  (13)  becomes 


(Au^  -  Au^*^)  [cosip^cosip^^  -  sinipgsintp^^] 


+  (Au  -  Au  ^)[sin®  cos®  ^  +  cos®  sinip 
X  X  o  y  o  y 


(16) 


+  (z  -  z,)[-cos®  sin<p  ^  -  sin®  cos<p  ^]Aq)  ^ 
b  o  y  o  y  y 


+  (x  -  X,  )  [-sin<p  sincp  ^  +  cos®  cosqj  '^jAtp  ^  =  0 
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(Au  -  Au  )cos(®  +  <p  )  +  (Au  -  Au  )sin(<p  +  w 
2z  oy  XX 


+  Acp^  ^(z  -  z,  )sin((|i  +  (P^  )  +  (x  -  X.  )cos(tp  +  <p^)/  =  0 
y  I  b  o  y  b  o  y  ; 


Next,  we  note  the  end  plane  does  not  rotate  about  the  local  Y-axis,  i.e.,  the 
will  be  no  component  of  rotation  about  vector  a^i*  A  small  rotation  about  the 
local  Y-axis  can  be  obtained  by  a  rotation  vector  with  respect  to  the  global 
system 


<Py.  <1'^) 

If  this  vector  is  normal  to  vector  a^ ,  then  it  has  no  rotational  component 
about  i.e. , 


£  *  li  =  *  [co8(<p  +  9,  ),  0, 

1  A  y  z  y  y 


•8in(®  +9  )]  =  0 

o  y 


9  cos(q)  +  9  -  9  sinC®  +  ®  ^)  =  0 

X  ’^o  y  z  ’^o  y 


The  incremental  form  of  Equation  (20)  becomes 


A®  cosC®  +  9  ^)  -  A<p  8in(9  +  9 

X  o  y  z  o  y 


^  I9  8in(9  +  9  +9  cos(9  +  9  ^)|  =  0 
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p  Finally,  the  essential  boundary  condition  is  that  the  rotation  q)y  of  all 

J  nodes  on  the  periphery  must  equal  the  rotation  about  the  local  Z-axis  ( a? 

[  vector  or  global  y-axis)  of  the  auxiliary  node,  i.e., 


I 


b 


and  in  incremental  form 


Affl  = 

y 


A<p 


(23) 


The  vertical  (global  x-axis)  rigid  body  motion  must  be  removed.  This  can  be 
accomplished  by  coupling  the  average  X-displacement  of  the  shell  nodes  on  the 
end  section  to  the  vertical  displacement  of  the  auxiliary  node.  This,  in  turn, 
is  set  to  zero.  For  element  i,  S8R  element  displacements  across  an  edge  are 
quadratic,  i.e.  (with  0  ^  5  ^  1)  the  vertical  displacement  u(5)  is  given  in 
terms  of  its  nodal  values  (u^^,  u^^-t-lt 


u(0  =  [(2^  -  1)(5  -  1),  4C(1  -  O,  C(2C  -  1)1 


k+1 


LV2J 


(24) 


For  uniform  finite  element  grid,  where  s  denotes  arc  length,  s£  arc  length  of 
element  i,  and  L  total  arc  length. 


£  . 
1 


(25) 


ds . 
1 


l^d^ 


(26) 


At  the  element  level  (element  i),  with  nodal  vertical  displacements  (u)^, 
'*k+l>  ^k+2) 


a 
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8" 

/ 


i 

u(t)d8 . 

X 


1 

u(c)l^d5 


\  ■ 

‘^k  +  l 

\+l 

V2 

(27) 


Since  there  are  N  elements  (with  2  corner  and  1  midside  nodes,  i.e.  3  nodes 
total),  there  will  be  2N-I-1  nodes  and  associated  nodal  values  of  vertical 
displacements.  Therefore,  by  Equation  (24) 


s=L 


/ 


u(s)ds 


i=N  5=1 


u(5)t.d5 

1 


6N 


[1, 


4, 


2.  4, 


{lx(2N+l) ) 


%+2N-l 

'^s+2N 

{(2N+1)  X  1) 


(28) 


Consequently,  the  boundary  condition 


L 


8=L 


/ 

s=0 


u(s)ds 


(29) 


can  be  recast , 

I 


in  view  of  Equation  (28),  to 


10 
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6N 


11. 


4, 


2.  4. 


2.  4,  1] 


8+1 


8+2N-1 


b 

u 

X 


(1  X  (2N+1)} 


|_“8+2N 


{(2N+1)  X  1) 


(30) 


for  a  uniform  grid,  where  the  nodal  vector  of  Equation  (30)  denotes  vertical 
displacements  of  shell  nodes  at  the  end  plane.  The  above  constraints  can  be 
incorporated  through  the  nonlinear  MPC  capability  of  ABAQUS. 1^2  Furthermore, 
from  Equations  (3)  and  (9),  Equation  (9)  can  be  rewritten  as 


u 

X 


b 

X 

o 


u  ^]8in(q)  +  9  ^) 

X  o  y 


(31) 


+ 


[z  + 
o 


u 


z 


U  ^]C08(<f  +  (p  )  *  0 

2  o  y 


Following  Reference  152,  we  eliminate  U2  from  Equation  (31)  first  (i.e.  at  the 
shell  node) 


u 

z 


b  bx^  ,  bx 

X  +  u  -  u  Hani?  +  «p  ; 
o  X  X  o  y 


Solving  for  <Px  from  Equation  (20) 


(32) 


(p  =  <p  tan(<p  +  <p  ^)  (33) 

X  z  o  y 


and  (Py  from  Equation  (22) 


<P 


y 


(34) 


Furthermore,  we  eliminate  the  vertical  translation  u^.  This  follows  from 
constraints  involving  the  elimination  of  u^.  All  previous  steps  are  explained 
in  Reference  152  or  can  be  found  in  the  FORTRAN  listing  for  MPC  constraints. 


»  XJk'i-d 


11 
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APPLICATIONS  AND  RESULTS 


Four  unstiffened  circular  cylinders  have  been  analyzed  using  the  finite 
element  program  ABAQUS^^^:  models  lOA,  16A,  20A  of  Reference  93,  and  20AI, 
which  is  an  imperfect  version  of  model  20A.  Note  that  they  are  comparatively 
thick  to  avoid  premature  collapse  and  to  exhibit  plasticity  effects  as  the 
applied  external  moment  increases  past  the  yield  moment.  The  above  models  fall 
in  the  range  of  "long  cylinders."  Their  perfect  discrete  analogues  would  fail 
by  ovalization  or  the  "Brazier  effect,"  or  by  plastic  deformation.  The 
experimentally  imperfect  models  failed  by  buckling  "failure"^^  as  displayed  on 
the  moment-curvature  curves. 

Table  1  gives  the  geometrical  and  material  properties  of  these  models  and 
the  stress-strain  curves  can  be  obtained  from  Reference  93.  In  Reference  93, 
the  external  loading  was  applied  by  the  use  of  a  so-called  "shear  span"  prior  to 
the  test  section  "bending  span."  To  create  the  external  moment  in  our  research, 
we  employed  an  additional  span  beyond  the  bending  span,  equal  to  the  "shear" 
span,  and  then  applied  a  fixed  angle  of  rotation.  Table  2  contains  some 
parameters  used  in  reducing  the  stresses,  moments,  curvatures,  and  forces  in 
nondimens ional  form,  and  they  can  be  used  to  determine  relative  magnitude  for 
critical  quantities  such  as  yield  moment,  etc. 

Mean  curvature  k  is  defined  as  the  ratio  of  the  sum  of  the  absolute  values 
of  direct  longitudinal  strains  at  0“  (top  of  cylinder,  i.e.  tension  side)  and 
180®  (bottom  of  cylinder,  i.e.  compression  side)  divided  by  the  undeformed 
diameter  of  the  shell.  Figures  3  through  5  are  the  relevant  M  vs.  k  curves  for 
models  lOA,  16A,  and  20A,  together  with  the  corresponding  experimental  results 
of  Reference  93.  Agreement  between  experimental  (dots)  and  computed  results  is 
good  for  all  three  cylinders.  Imperfection  sensitivity  and  the  analysis  of 
medium  and  short  cylinders,  where  short-wave  length  buckling  may  control  the 
collapse  mechanism,  will  be  addressed  elsewhere. 

Figures  6  through  8  display  the  work  done  versus  the  moment  parameter  u  = 
Mext/^0*  where  Mq  =  4R^  hOy,  R  =  mean  radius  of  cylinder,  E  =  Young's  Modulus, 
Mext  ~  external  bending  moment,  h  =  shell  thickness,  and  Oy  =  material  yield 
stress . 

We  define  two  parameters,  longitudinal  stress/stress  at  bifurcation  and 
hoop  stress/stress  at  bifurcation,  where  stress  at  bifurcation  is 


o 


CR 


h 

R 


Figures  9  through  11  represent  plots  of  the  longitudinal  inplane  stress 
parameter  versus  angular  position.  Additional  local  bending  stresses  across  the 
shell  thickness  are  not  addressed.  Note,  however,  that  the  longitudinal 
membrane  stress  distribution  in  Figures  9  through  11  does  not  follow  simple  beam 
theory.  The  stress  distribution  from  the  0°  and  18U“  points  (top  and  bottom  of 


TABLE  1.  GEOMETRICAL  AND  MATERIAL  PROPERTIES  OF  MODELS 


BENDING  MOMENT  -  MEAN  CURVATURE 
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BENDING  MOMENT  -  MEAN  CURVATURE 


(X10~3) 

I  I  I _ 1  i _ 1 _ L 

0.8  1.6  2.4  3.2  4.0  4.8  5.6 

MEAN  CURVATURE  (1/FT) 


FIGURE  5.  MOMENT  -  CURVATURE  PLOTS  FOR  MODEL  20A 
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EXTERNAL  WORK  -  BENDING  MOMENT/M, jljiivi  MOMENT 


20  40  60  80  100  120  140  160 


WORK  DONE  (KIP-lN) 

FIGURE  6.  M/Myijliy,  WORK  DONE  (KIP-IN)  FOR  MODEL  lOA 
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FIGURE  10.  LONGITUDINAL  STRESS  (STEP  31)  VERSUS  ANGULAR  POSITION  FOR  MODEL  16A 
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the  halt-section)  is  reduced  (.tension  zone)  or  increased  (compression  zone) 
slightly  before  it  assumes  a  linear  form. 

Figures  12  through  14  display  the  variation  of  the  hoop  stress  parameter 
with  angular  location.  The  distribution  of  this  stress  over  the  half-section 
clearly  indicates  maximum  compressive  hoop  stresses  around  the  90“  location, 
with  corresponding  maximum  tensile  hoop  stresses  at  0“  and  180“,  respectively. 
These  stresses,  however,  are  smaller  than  the  longitudinal  membrane  stresses  by 
an  order  of  magnitude.  All  stresses  shown  correspond  to  the  final  point  of  the 
moment-curvature  plots. 

Close  examination  of  Figures  15  through  17  reveals  that  the  vs.  R/Rcr 

curves  have  a  slope  at  the  origin  of  approximately  1.0.  This  agrees  fairly  well 
witn  the  initial  slope  predicted  by  Reissner's  nonlinear  theory^^  as  well  as 
Von  Karman ' s  linear  analysis.^  In  the  present  notation,  Reissner's  relation¬ 
ship  between  M/Mcr  and  k/k(;^g  can  be  written  as 

M/Mcr  =  k/k^j^[1.0  -  0.5  (k/k^j^)^  -  (1/6)  (k/k^^^)^  .  .  .  .  ]  (35) 


with  an  obvious  slope  of  1.0  at  the  origin,  maximum  value  of  M/M(jg  -  0.5011 
at  k/kcR  =  0.7  1954.  The  maximum  value  of  0.5011  is  much  higher  than  both 
experimental  and  numerical  results,  showing  the  effects  of  plasticity  which 
Reissner's  theory  does  not  account  for.  Note  that  plots  of  p  =  (MgxT^^O^ 
versus  work  done  (Figures  6  through  8)  clearly  display  that  ultimate  moment  is 
reached  at  p  =  1  as  predicted  by  beam  plastic  theory.  Actually,  it  exceeds 
p  =  1  by  about  10  percent. 

Figure  18  displays  the  deformed  and  undeformed  profiles  of  cylindrical 
model  20A.  Figure  19  represents  the  corresponding  deformed  and  undeformed 
profiles  of  the  imperfect  version  of  model  20A,  which  from  now  on  will  be 
referred  to  as  model  20AI.  Figure  20  is  an  enlargement  of  the  cross-sections  at 
uiidlength  before  and  after  deformation  of  model  20AI.  Figure  21  shows  the 
moment-curvature  plots  of  model  20A1  compared  with  the  experimental  data 
(circled  points).  Similarly,  Figures  22  through  25  show  results  for  model  20AI 
of  the  following  pairs  of  variables: 


P  =  Mrxt/^O  ^/^ultIMATE^  ~  Work  Done  (kip-in.) 


LONG  , 

0  /o 

INPLANE  CR 


Angular  Position 


liOOP^^  Angular  Position  and  M/M(;r  ~  k/k(jR 
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FIGURE  18.  CYLINDRICAL  SHELL  (MODEL  20A)  SUBJECT  TO  END  BENDING  MOMENT  VIEWED  FROM 
100  MOO ',500  "  AT  STEP  3,  INCREMENT  55 
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ri  INC  £Fi  WITH  IMPcRF  tl‘  ^  I  I.'N:.  ■  h.\tl  J  TO  ENj  rClIJPLE 

«  ittfacHCBT  2 


figure  19.  IMPERFECT  CYLINDRICAL  SHELL  (VARIANT  OF  MODEL  20A)  SUBJECT  TO  END  BENDING 
MOMENT  VIEWED  FROM  100 MOO  '.SOO  '  AT  STEP  6,  INCREMENT  2 


VER 


FIGURE  20.  UNDEFORMED  AND  DEFORMED  HALF  CROSS  SECTION  OF  IMPERFECT  VERSION  OF 
MODEL  20A  AT  MIDLENGTH 


BENDING  MOMENT  -  MEAN  CURVATURE 


MEAN  CURVATURE  (1/FT) 


EXTERNAL  WORK  -  BENDING  MOMENT/ULTIM  MOMENT 


FIGURE  22.  M/MULTIM  --  WORK  DONE  <KIP-IN.)  FOR  IMPERFECT  VERSION  OF  MODEL  20A 
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FIGURE  24.  HOOP  STRESS  DISTRIBUTION  (STEP  6,  INCREMENT  2,  WHICH  CORRESPONDS  TO  FINAL 
PLOTTED  POINT  ON  FIGURE  21)  VERSUS  ANGULAR  POSITION  FOR  IMPERFECT  VERSION 
OF  MODEL  20A 
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At  this  point,  note  that  the  relevant  figures  corresponding  to  the  perfect  model 
20A  are  8,  11,  14,  and  17,  respectively.  The  imperfection  in  the  radial 
displacement  of  the  original  surface  from  the  mean  radius  R  of  imperfect  model 
20AI  varied  according  to  the  formula 


h  sin  cos(  10  0) 

Imperfect  model  20AI  was  generated  with  one-half  wave  axially  and  10  waves 
peripherally . 

Figure  26  presents  superimposed  moment-curvature  curves  for  both  models  20A 
and  20A1.  Notice  that  since  they  are  comparatively  thick  and  fail  by 
plastif ication,  unlike  ovalization  or  bifurcation  failure,  they  are  not 
imperfection  sensitive.  The  response  of  the  imperfect  model  20AI  is  very 
similar  to  the  perfect  one. 

Table  3  summarizes  digitized  results  pertaining  to  bending  moment,  angle  of 
rotation  at  the  end,  where  the  external  load  is  applied,  as  well  as  strain 
energy,  work  done,  and  plastic  dissipation  for  model  20A.  Note  that  because  of 
an  existing  error  in  ABAQUS,  concerning  how  energies  are  computed,  strain  energy 
and  plastic  dissipation  do  not  agree  exactly  with  the  work  done  (0.120489  + 
0.290966  =  0.411455  compared  with  0.417864).  In  fact,  this  difference  decreases 
as  we  march  along  the  load-deformation  curve.  Table  4  gives  the  stress 
distribution  for  load  step  20  of  model  20A.  Tables  5  and  6  give  the 
corresponding  information  for  model  20AI,  and  Tables  7  and  8  give  the  actual 
points  plotted  in  Figure  26  for  both  models  20A  and  20A1 . 

Finally,  a  fine  and  extremely  important  point  pertaining  to  the  modeling 
issue  must  be  addressed.  The  experimental  set  up  involved  the  analysis  of  a 
three  span  beam  (in  all  cases)  subject  to  vertical  self-equilibrating  shear 
loads,  to  simulate  overall  bending  over  the  two  end  spans.  The  middle  span  had 
no  loads  applied.  In  this  analysis ,  one  end  span  and  half  of  the  center  span 
(symmetry)  were  modeled.  A  rotation  was  enforced  through  the  "auxiliary  node" 
concept  and  the  MFC  constraints.  Subsequent  computations  without  the 
"additional"  end  span  gave  a  slightly  different  response. 

Consequently,  what  constitutes  an  adequate  additional  span  to  be  included 
for  proper  response  has  not  been  determined  but  is  of  major  importance  in  this 
analysis.  In  addition,  to  complete  studies  on  the  bend-buckling  modeling  of 
cylindrical  shells,  "short"  and  "medium"  length  tubes  with  or  without  ring 
stiffeners  must  be  addressed. 

In  closing  this  discussion  it  must  be  stressed  that,  by  neglecting  inertia 
effects,  a  reasonably  successful  first  approximation  of  the  problem  of  a  long 
straight  circular  cylinder  shell  subject  to  end  couples  has  been  developed. 


NSWC  TR  84-380 


SEND  INC  MOMENT  VERSUS  REAM  CLJflVATLIRE 


FIGURE  26.  BENDING  MOMENT  DISTRIBUTION  VERSUS  MEAN  CURVATURE  FOR  BOTH  PERFECT 
AND  IMPERFECT  VERSIONS  OF  MODEL  20A 


TABLE  3.  POST-PROCESSING  INFORMATION  FROM  ABACUS  FOR  MODEL  20A 


TABLE  4.  POST-PROCESSING  INFORMATION  (STRESS  DISTRIBUTION)  FROM  ABACUS  FOR  MODEL  20A 
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TABLE  6.  POST- PROCESSING  INFORMATION  (STRESS  DISTRIBUTION)  FROM  ABACUS  FOR  IMPERFECT  VERSION 
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MOMENT-CURVATURE  RESULTS  FOR  PERFEC^  VERSION  OF  MODEL  20A 
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TABLE  8.  MOMENT-CURVATURE  RESULTS  FOR  IMPERFECT  VERSION  OF  MODEL  20A  (MODEL  20AI) 


CURVATURE  K 

MOMENT  M 

(1/FT) 

(KIP-FT) 
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0. 3379i£93E~03 

0, 52431121E+02 

0.  l£75  4.641E-0£ 

0,  197A49i.0E  +  03 

0.  -'cilciYioUE-OE 

0.  4  3 /  Y  Y0i  lE-i'03 

0 , 3190451 3E-0E 

0. 59942395 E+03 

V  5  3  i ^  1  5  El  ^  ^ 

0. 405  ■^453  3E  +  03 

0. 3439 35 355-05 

0. 410^S504E+O3 

5. 56557399E-05 

0 . 4155 135EE+03 
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0. 4195935 lE+03 
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0.  495  325  5  9E--02 
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SUMMARY 


A  modeling  strategy  is  established  to  obtain  the  moment-curvature  relation 
as  well  as  the  relation  between  the  work  done  and  the  applied  moment  for 
circular  cylindrical  shells.  This  is  achieved  by  using  the  nonlinear  finite 
element  program  ABAQUS  in  conjunction  with  preprocessing  and  postprocessing 
computer  programs.  The  results  compare  favorably  with  experimental  curves 
reported  in  the  open  literature.  Such  analysis  is  of  potential  use  in 
predicting  critical  bending  moments,  ultimate  moment  for  ship  hulls,  pipe  bends 
in  nuclear  reactors,  submarine  pipelines,  etc. 
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NOMENCLATURE 


11  >12  .13 


E 

FO 

l^R 

h 

I 

k 


L 


H 

Mbr 

Mcr 

Mo 

N 


=  Local  vectors  along  the  local  Y,  Z,  and  X  axes.  Local  z-axis 
is  parallel  to  global  y.  Initially,  local  y  is  not  parallel, 
but  is  coplanar  with  global  x  and  in  the  same  direction. 
Initially,  local  X  is  coplanar  with  global  z  and  points  in  the 
same  direction  (Figure  2). 

=  Young's  Modulus 

=  Axial  force  based  on  yield  stress 

=  Critical  curvature  at  bifurcation  (see  Table  2) 

=  Shell  thickness 

=  Moment  of  inertia  of  undeformed  section  (Table  2) 

=  Subscript  in  nodal  values  uj^,  uj^+j^,  ui^+2  (Eq.  (24)). 

Mean  curvature  defined  as  the  sum  of  absolute  value  of  direct 
axial  strains  at  the  top  and  bottom  of  the  cylinder  divided  by 
the  undeformed  diameter  of  the  shell. 

=  Total  length  of  cylindrical  shell  (Figure  1).  In  Equation 

(25),  L  is  total  arc  length  of  end  section  ABC  (Figure  2)  prior 
to  deformation. 

=  Length  of  side  of  element  i  on  arc  where  end  rotation  is 
applied  (undeformed  side  ABC  on  Figure  2) 

=  Critical  bending  moment  by  Brazier  (Table  2) 

=  Critical  bending  moment  of  bifurcation  (Table  2) 

=  External  bending  moment  calculated  in  the  form  of  a  reaction  in 
this  analysis 

=  Plastic  moment  or  ultimate  moment  based  on  yield  stress 
(Table  2) 

=  Total  number  of  equal  elements  on  undeformed  arc  ABC  (Figure  2) 
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NOMENCLATURE.  (Cont.  ) 


=  Final  position  vector  on  arc  A'B'c'  (Figure  2,  after 
deformation) 


r 


o 


-  Initial  position  vector  of  any  point  on  periphery  of  shell  (arc 
ABC  on  Figure  2)  at  which  the  end  rotation  will  be  applied 


rb 

-Lo 


=  Initial  position  vector  of  auxiliary  node 


R 


=  Mean  shell  radius 


u  =  In  the  "APPLIED  LOADING"  section,  u  is  the  vector  displacement 

after  deformation;  it  is  not  to  be  confused  with  u 

ub  =  In  the  context  of  analysis  given  in  "APPLIED  LOADING,"  vector 

displacement  of  auxiliary  node  after  deformation;  it  is  not  to 
be  confused  with  u 


Ux,Uy,U2  =  Vertical,  transverse  and  longitudinal  translations  (along 

global  X,  Y,  Z  axes).  Also  referred  to  as  u,  v,  and  w. 


(u) 

Vl>V2 


u 


s 


u 

V 

w 


V 


=  Vertical  displacement  at  any  point  ^  on  an  element  side 

=  Nodal  values  of  displacement  u(£))  at  end  node  k,  midside  node 
k+1,  end  node  k+2. 

=  Nodal  value  at  s  (point  on  arc  where  end  rotation  is  being 
applied) 

=  Translation  along  global  x-axis 

=  Translation  along  global  y-axis 

=  Translation  along  global  z-axis 

=  Angle  from  vertical  global  x-axis  up  to  end  of  arc  to  be 
analyzed 

=  Poisson's  ratio 


K  =  Local  normalized  variable  in  interval  (0,  1)  employed  in 

describing  vertical  displacement  distribution  u(^). 

=  Critical  stress  at  bifurcation  (Table  2) 

oy  =  Yield  stress  of  material 

<pQ  =  Initial  inclination  of  cylindrical  shell  about  y-axis.  For 

straight  tubes  (present  case)  (pQ  =  0 

=  Rotation  about  global  x-axis 
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